H II regions: Witnesses to massive star formation 
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ABSTRACT 



We describe the first three-dimensional simulation of the gravitational collapse 
of a massive, rotating molecular cloud that includes heating by both non-ionizing 
and ionizing radiation. These models were performed with the FLASH code, 
incorporating a hybrid, long characteristic, ray tracing technique. We find that 
as the first protostars gain sufficient mass to ionize the accretion fiow, their H ll 
regions are initially gravitationally trapped, but soon begin to rapidly fluctuate 
between trapped and extended states, in agreement with observations. Over 
time, the same ultracompact H ll region can expand anisotropically, contract 
again, and take on any of the observed morphological classes. In their extended 
phases, expanding H II regions drive bipolar neutral outflows characteristic of 
high-mass star formation. The total lifetime of H II regions is given by the global 
accretion timescale, rather than their short internal sound-crossing time. This 
explains the observed number statistics. The pressure of the hot, ionized gas does 
not terminate accretion. Instead the flnal stellar mass is set by fragmentation- 
induced starvation. Local gravitational instabilities in the accretion flow lead to 
the build-up of a small cluster of stars, all with relatively high masses due to 
heating from accretion radiation. These companions subsequently compete with 
the initial high-mass star for the same common gas reservoir and limit its mass 
growth. This is contrary to the classical competitive accretion model, where the 
massive stars are never hindered in growth by the low-mass stars in the cluster. 
Our flndings show that the most signiflcant differences between the formation of 
low-mass and high-mass stars are all explained as the result of rapid accretion 
within a dense, gravitationally unstable, ionized flow. 
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Introduction 



Massive stars influence the surrounding universe far out of proportion to their numbers 
through ionizing radiation, supernova explosions, and heavy element production. Their 
formation requires the collapse of massive interstellar gas clouds with accretion rates 
exceeding 10~^ yr~^ ^ Beuther et al.|[2b02 ; Belt ran et al.|[2006 ) to reach their flnal masses 



before exhausting their nuclear fuel (Keto & Wood 2006). Massive stars can ionize the 



gas around them, forming H II regions, that traditionally have been modelled as spherical 



bubbles expanding in a uniform medium (Spitzer 1978). This approach, however, fails to 



explain their observed numbers and morphologies (Wood & Churchwell 1989; Kurtz et al. 



1994). A more recent alternative picture models H ll regions as the ionized section of the 



accretion flow that feeds the young massive star (Keto 2002, 2007). The rich diversity of 
observed H ll regions thus bears witness to the complexity of high-mass star formation. 

H II regions form around accreting protostars once they exceed ^ 10 M©, equivalent 
to a spectral type of early B. Thus, accretion and ionization must occur together in the 
formation of massive stars. The pressure of the 10^ K ionized gas far exceeds that in the 



10 K accreting molecular gas, producing feedback in the form of ionized outflows (Keto 



2002, 2003, 2007). 



High-mass stars form in denser and more massive cloud cores (Motte et al. 2008) 



than their low-mass counterparts (Myers et al. 1986). High densities also result in local 



gravitational instabilities in the accretion flow, resulting in the formation of multiple 



additional stars (Klessen & Burkert 2000; Kratter & Matzner 2006). Young massive stars 



are almost always observed to have companions (Ho & Haschick 1981), and the number 



of their companions signiflcantly exceeds those of low-mass stars (Zinnecker & Yorke 



2007). Such companions influence subsequent accretion onto the initial star (Krumholz 



et al. 2009a). Observations show an upper mass limit of about 100 M©. It remains 
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unclear whether hmits on internal stability or termination of accretion by stellar feedback 



determines the value of the upper mass limit (Zinnecker & Yorke 2007). 



Around the most luminous stars, the outward radiation pressure can counterbalance the 
inward gravitational attraction. A spherically symmetric calculation of radiation pressure 



on dust yields equality at just under 10 Mq (Wolfire & Cassinelli 1987). However, the 
dust opacity is wavelength dependent, the accretion is non-spherical, the mass-luminosity 
ratio is different for multiple companions than for a single star, and the momentum of the 



accretion flow must be reversed, rather than just achieving static force balance (Larson & 



Starrfield 


1971; 


Kahn 


1974; 


Yorke & Kriigel 


1977; 


Nakano et al. 


1995 



Observations (Chini et al. 2004; Patel et al. 2005; Keto & Wood 2006; Beltran et al. 2006) 



provide evidence for the presence of all these mitigating factors, and numerical experiments 



combining some of these effects (Yorke & Sonnhalter 2002; Krumholz et al. 2007b) conflrm 



their effectiveness, showing that radiation pressure is not dynamically signiflcant below the 
Eddington limit. 

The most significant differences between massive star formation and low-mass star 
formation seem to be the clustered nature of star formation in dense accretion flows and the 
ionization of these flows. We present the flrst three-dimensional simulations of the collapse 
of a molecular cloud to form a cluster of massive stars that include ionization feedback, 
allowing us to study these effects simultaneously. During recent years the problem of 
massive star formation was tackled in a number of three-dimensional numerical studies (see 
e.g. 



Klessen et al. 2009, for a recent review). However, none of these simulations included 



the self-consistent ionization feedback from massive protostars that dynamically condense 



out of gravitationally unstable regions. Dale et al. (2005) performed the flrst hydrodynamic 



calculation of a star cluster forming region to include a simplifled treatment of ionization 
feeback, but only from a single source that was inserted into a pre-existing simulation 



- 5 - 



and whose ionizing flux did not depend on the mass of the growing (proto-)star. Also, 
heating by the accretion luminosity and from the stellar photosphere is neglected in these 
calculations, but is included in ours. Other calculations that incorporate ionizing radiation 
focus on triggered star formation and turbulence by external sources ( [Dale et al.||2007a|[^ 



Gritschneder et al. 2009). Using a flux-limited diffusion approximation, Krumholz et al 



(2007b) incorporated radiation feedback in their calculations of massive star formation. 



These calculations showed that radiation pressure from non-ionizing radiation cannot halt 
accretion onto the massive protostar and therefore is not able to set an upper mass limit 



of such stars (see also Krumholz et al. 2009a). Nevertheless, these simulations lack the 
important feedback from ionizing radiation distinctive of massive star formation. 



In the following, we present our numerical method (§ [2]) and give a detailed analysis 
of the simulation results. We discuss the accretion history of the stellar cluster (§ [3]), the 
fragmentation of the rotationally flattened structure, and the generation of ionization-driven 
bipolar outflows (§|4]), and compare synthetic radio continuum and line observations with 
observed data (§ [s]). Finally, we discuss the relation of our numerical model to previous 
work (§ IgI), and the relevance of our flndings for massive star formation (§ [T]). 



2. Numerical Method 



2.1. Algorithm 



We use a modifled version of the FLASH code (Fryxell et al. 2000), an adaptive-mesh 



code that integrates the compressible gas dynamic equations with self-gravity and radiation 
feedback. Our modiflcations include the propagation of both ionizing and non-ionizing 
radiation to follow heating, though not radiation pressure, a reflnement criterion to 



resolve the Jeans length ( [Banerjee et al.||2004p , various cooling processes of relevance for 
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protostellar collapse (Banerjee et al. 2004; Banerjee & Pudritz 2006; Banerjee et al. 2006); 



and Lagrangian sink particles (Banerjee et al. 2009; Federrath et al. 2009). Here we report 



on our new modifications to the code. 



2.1.1. Radiation 

To propagate radiation, we use the hybrid characteristics raytracing module originally 
developed by Rijkhorst et al.| ( |2006 ). Parallel raytracing on a block-structured adaptive 



mesh requires the solution of an inverse problem: Given the position in the domain, the 
corresponding block identifier is needed. This problem was originally approached by storing 
the block identifier in one single large array corresponding to a fully refined domain. Each 
entry in this array could be mapped to a (potential) cell in the domain and vice versa, so 
that the block identifier could be obtained by looking up the stored value. This solution 
is straightforward and fast, but it violates the principles of adaptive mesh refinement 
and becomes extremely memory consuming at high resolution. Collapse simulations are 
impossible using this method. 

Instead of creating one large array, we use the fact that FLASH stores hierarchical 
information about its block structure. Every block in the adaptive-mesh hierarchy has 
information about its parent and child blocks as well as on its neighbors at the same level of 
refinement. Valid data is stored for blocks of highest refinement only, which are called leaf 
blocks. Neighboring leaf blocks can diflFer only by one level of refinement. This hierarchical 
information makes it possible to perform a tree walk in the adaptive-mesh hierarchy. 

Starting from a block whose identifier is known, we determine the direction in which 
the ray leaves the block. If there is a neighboring block on the same level of refinement, 
it can either still be a leaf block or have one level of children. In the former case the new 
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block is already found, in the latter case the child blocks must be checked. If there is no 
neighboring block at the same refinement level, the new block must be the neighbor of the 
parent of the current block in this direction. 

Although the method is more complicated than the original one, it is equally fast since 
the hierarchical data can be eflFectively stored and communicated among processors. 

The radiation module has been tested for its capability of accurately casting shadows 
(Rijkhorst et al. 2006)) and tracing R-type ionization fronts in a cosmological setting (Iliev 



et al. 2006). For simulations of massive star formation, the expansion of D-type ionization 



fronts should be adequately modeled. We have tested the ionization physics in our model 



against an approximate solution found by Spitzer (1978) for the expansion of a D-type 



ionization front into a homogeneous medium. Given the Stromgren radius Rs and the 
sound speed of the ionized gas Cg, the radius R of the ionization front at time t is given by 

A comparison of our numerical simulation with this analytic solution is plotted in 
Fig. [ij We show the data from two different simulations, both of which had initial density 
p = 3.3 X 10~^^ g cm~^, ionizing luminosity = 1.4 x 10^^ s~^, and were run with a 
maximum resolution of 0.25 pc. In the first run, we used an isothermal equation of state 
(7 = 1) and included no cooling processes. In the second, we set 7 = 5/3 and used the 



cooling curve from Dalgarno & McCray (1972). In both simulations, the ionization fraction 



and temperature of the ionized gas was determined self-consistently by solving a rate 



equation for the ionization fraction as discussed in Rijkhorst et al. (2006). 



We compare the numerical result in each case with the result from equation ([T]) using 
the appropriate sound speed in the ionized gas Cg = {^ykTi/ /llY^'^ ^ where k is Boltzmann's 
constant and /x is the mean mass per particle. The agreement of each run with the analytical 



-8- 



solution and with each other is acceptable. The decent agreement even in the presence of 
cooling demonstrates that numerical diffusion from the cold, dense, neutral shell into the 
hot, ionized interior does not lead to significant unphysical cooling of the hot interior. 

We show a slice of the expanding H ll region in figure [2} The whole box shown has 
an effective numerical resolution of 128^ grid points. Small deviations from a perfectly 
spherical shape are visible. This is because the sphere is initially not very well resolved, 
so that dynamical processes inside the H ll region can lead to an amplification of grid 
effects. Resolution studies in two dimensions that we have performed show that this effect 
diminishes with increasing resolution. 



2.1.2. Protostellar model 
To trace the collapse of individual fragments, we use a sink particle technique that 



we developed for the FLASH code (Federrath et al. 2009) (see also Bate et al. (1995)). 



Here, sink particles are created if the local density exceeds a critical value pcrit and the 
region within the sink particle radius rgmk is gravitationally bound and collapsing. For 
the high resolution simulations (Run A and Run B), we use pcrit = 7 x 10~-^^gcm~^ 
and Tgink = 590 AU, for the lower resolution simulations (Run Ca and Run Cb), we take 
Pcrit = 1 X 10~^^gcm~^ and rgmk = 1319 AU. The sink particles gain mass through accretion 
of overdense gas within the accretion radius that is gravitationally bound to it. The 
particles interact gravitationally with the gas and with themselves and are free to move 
within the simulation box, independent of the underlying grid structure, i.e. they are 
treated fully Lagrangian. 

We use the properties of these Lagrangian sink particles as sources for the radiation 
module. A prestellar model is used to determine luminosity and effective temperature of 
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12 3 4 

t (Myr) 

Fig. 1. — Expansion of a D-type ionization front in a homogeneous medium. The plots shows 
the radius R of the ionization front as function of time t. The analytical data is dashed, the 
numerical data is solid. The plot shows a run with an isothermal equation of state (7 = 1) 
without cooling and a run with 7 = 5/3 and a cooling curve. The results agree well, both 
with the analytic solution, and with each other. 



^ 10 ^ 




Fig. 2. — Log density slice of an ionization front in a homogeneous medium. The sheU devi- 
ates a httle from a perfect sphere due to waves running through the ionized gas, amphfying 
grid effects. The deviations from the spherical shape decrease with increasing resolution. 
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the source as function of protostellar mass and accretion rate. We set the steUar luminosity 



by a zero- age main sequence (Paxton 2004j) and the accretion luminosity by interpolation 



from a more detailed model (Hosokawa & Omukai 2009). 



The photoionization rate and the photoionization heating rate are set by the specific 
mean intensity along a ray 



^ /; rj. — X -exp-Tion r . 

V r / 2c2 exp{hiy/kBTstaT) - 1 

Here, r is the distance from the source, Tgtar the radius of the star, Tgtar the effective 
temperature of the star, u the frequency, c the speed of light, h Planck's constant, /cb 
Boltzmann's constant and rion the optical depth for the ionizing radiation, which is 



(2) 



calculated using the absorption coefficient of atomic hydrogen (Rijkhorst et al. 2006). 



Hence, the strength of the ionizing radiation is totally determined by Tgtar and Tgtar- We use 
a zero-age main sequence (ZAMS) model derived from the freely available stellar evolution 



code EZ (Paxton 2004) to calculate these quantities. 



To also include heating by non-ionizing radiation, we have extended the ray tracing 



module to calculate Planck mean opacities following the Krumholz et al. (2007b) 



interpolation of data from [Pollack et aL (1994). We use the non-relativistic approximation 



for the dust heating term (see e.g., Krumholz et al. 2007a) 



tvppcu 



(3) 



with the Planck mean opacity /^p, the gas density p and the total radiation energy u. In the 
raytracing approximation, the heating due to stellar radiation is thus given by 



St I 



star 



(4) 



with the Stefan-Boltzmann constant a and the Planck mean opacity rp^ 



In addition to the stellar heating, we also include a model for accretion heating. To 
this end, we assume that the potential energy of the gas which is released during accretion 
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is fully converted into radiation at the accretion radius race, so that the accretion luminosity 
is given by 

^acc G (5) 

^acc 

with Newton's constant G, the protostellar mass M and the accretion rate M . Since we 
do not account for diffuse radiation, the accretion heating also originates from the sink 
particles and leads to another heating term of the form 

racc(r) = a (^)^p[T(r)]p(r) exp[-rp^(r)]r,te. (6) 

To determine the accretion radius race and the effective temperature T^^c of the accretion 
luminosity, we interpolate the results of the prestellar evolution model of |Hosokawa fc 



Omukai (2009) from the current accretion rate and protostellar mass. 



2.2. Initial Conditions 

We start our calculations with a molecular cloud with a mass of 1000 Mq and an initial 
temperature of 30 K. The cloud has a flat inner region with 0.5 pc radius, surrounded by a 
region in which the density drops as r~^/^. The core density is 1.27 x 10~^^ gcm~^. The 
whole simulation box has an edge length of 3.89 pc. The core is in solid body rotation with 
a ratio of rotational to gravitational energy (5 = 0.05. We then run the simulation allowing 
sink particles to continue forming for 130 kyr after the formation of the first one. Note that 
the total simulation time is 0.75 Myr. 

These initial conditions represent an extrapolation from lower mass prestellar cores, 
which are well described by density profiles similar to Bonnor-Ebert spheres (see e.g.. 



Pirogov 2009), towards higher masses and larger scales. After a simulation runtime of 
500 kyr, the flat inner part has totally disappeared, and a centrally condensed structure 
with a peak density of 10~^^ g cm~^ and a radius of 0.65 pc has emerged. This is the stage 
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at which infrared dark clouds are typicaUy observed (see e.g., Beuther & Henning 2009). 



In some runs, we choose to suppress secondary sink formation and instead use a 
density-dependent temperature floor to prevent runaway coUapse of dense blobs of gas. We 



need to resolve the local Jeans length with n > 4 cells (Truelove et al. 1997). To do so, we 
introduce a dynamical temperature floor 



T 

J. rr 



-p{nAxy 



(7) 



where fi is the mean molecular weight, the proton mass, and Ax the cell size. 



We consider two full resolution simulations with a cell size at the highest reflnement 
level of 98 AU (11 reflnement levels). In one. Run A, we artiflcially suppress disk 
fragmentation using the dynamical temperature floor, while in the other. Run B, we permit 
it and dynamically form secondary sink particles. In addition, we ran two half resolution 
simulations with maximum resolution of 196 AU (10 reflnement levels) and suppressed 
secondary sink formation. One of these runs. Run Ca, starts with identical initial conditions 
to Runs A and B, while Run Cb has an additional m = 2-perturbation of 10% of the gas 



density (Boss & Bodenheimer 1979). The high-resolution stellar cluster simulation Run B 
reaches ^ 7.2 x 10^ grid cells by the end of the simulation, and required 2.5 x 10^ CPU 
hours to complete. An overview of the runs with their diflFerent model parameters is given 
in Tab. [2 



Name 


Resolution 


Multiple Sinks 


Run A 


98 AU 


no 


Run B 


98 AU 


yes 


Run Ca 


196 AU 


no 


Run Cb 


196 AU 


no 



Table 1: Overview of collapse simulations. 
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3. Accretion 

We first examine tlie accretion liistories of our different models. Figure [3] shows that in 
Run A, with only one sink particle allowed to form, nothing halts accretion onto the central 
protostar. It continues to grow at an average rate of M ^ 5.9 x lO~^M0yr~^ until we stop 
the calculation when the star has reached 72 Mq. The increasingly massive star ionizes the 
surrounding gas, raising it to high pressure. This gas breaks out above and below the disk 
plane, but it cannot halt mass growth through the disk mid-plane . 

We ran our lower resolution simulations with a single sink particle. Runs Ca and Cb, 
until the sink particles reached masses of 94 and 100 Mq without any reduction in the 
mass accretion rate. Figure [4] shows the accretion history for these simulations. One can 
clearly see that neither the additional perturbation nor the change in resolution has any 
significant influence on the overall accretion behavior. Runs Ca and Cb accrete at mean 
rates of 4.6 x 10~^ Mq yr~^ and 5.8 x 10~^ Mq yr~^, almost identical to the value for Run A, 
showing that this quantity is well converged. These simulations with a single sink particle 
strongly suggest that ionization feedback alone does not limit massive star accretion. 

In contrast, in Run B, where multiple sink particles are allowed to form, two subsequent 
sink particles form and begin accreting soon after the first one, and many more follow 
within the next 10^ yr (see Fig. [s]). By that time the first sink has accreted 8Mq. Within 
another 3 x 10^ yr seven further fragments have formed, with masses ranging from 0.3 
to 4.4 M0 while the first three sink particles have masses between 10 Mq and 20 M©, all 
within a radius of 0.1 pc from the most massive object. Accretion by the secondary sinks 
terminates the mass growth of the central massive sinks. Material that moves inwards 
through the disk driven by gravitational torques accretes preferentially onto the sinks at 



larger radii (Bate 2002). Eventually, hardly any gas makes it all the way to the center to 



fall onto the most massive objects. This fragmentation-induced starvation prevents any star 



1000 



10-2 



Run A 
Run B 

Run B (sum) 



Run A 
Run B (1^*) 




iL 
0.60 



^10-^ 




0.65 



0.70 



0.75 



0.60 



0.65 0.70 



0.75 



t(Myr) 



t(Myr) 



Fig. 3. — Accretion history of the single (Run A) and multiple (Run B) sink simulations. 
Run A was stopped at 72 Mq, while no sink particle in Run B exceeds 25 Mq over the 
simulation runtime. The left hand plot shows the sink particle masses for Run A (green), 
the individual sink masses of Run B (blue) as well as the total mass in sink particles in Run 
B (red) . The right hand plot shows the accretion rates of the sink particle in Run A (green) 
and sink particle which forms first in Run B (blue), which also turns out to end up as the 
most massive at the end of the simulation. The most massive stars in Run B are those which 
form early and keep accreting at a high accretion rate. While the accretion rate in Run A 
never drops below 10~^ yr~^, accretion onto the most massive sink can drop significantly 
below this value and even be stopped totally in Run B. 
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0.60 0.65 0.70 0.75 0.60 0.65 0.70 0.75 

t(Myr) ^(Myr) 



Fig. 4. — Accretion history of the high resolution (Run A) and low resolution simulations 
with the same initial conditions (Run Ca), and with an additional initial perturbation (Run 
Cb). The accretion process proceeds similarly in all three simulations. Run A was stopped 
at 72 M0, Run Ca at 100 Mq and Run Cb at 94 M©. The mean accretion rates of 5.9 x 
10""^ M0 yr~^, 4.6 x 10~^ yr~^ and 5.8 x 10~^ Mq yr~^, respectively, are also very similar. 
There is no indication that ionization may stop accretion. 
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from reaching a mass > 25 Mq in this case. 

The accretion behavior in Run B contrasts sharply with competitive accretion models 



(Bonnell et al. 2001, 2004), which have no mechanism to turn off accretion onto the most 
massive stars. In these models, material falls all the way down to the massive stars sitting 
in the center of the gravitational potential which thereby take away the gas from the 
surrounding low-mass stars. In our fragmentation-induced starvation scenario, exactly the 
opposite happens. 

The different behaviour of the multiple sink simulation (Run B) and the single sink 
simulation (Run A) can be understood by looking at slices of density in the disk plane, as 
shown in Figures [5] and (H) In Run B, secondary sink particles form in the disk filaments, 
which remove gas from the common reservoir. Consequently, the most massive stars cannot 
maintain their high accretion rates. It is this halting of the accretion flow that allows the 
ionizing radiation to create a bubble around the massive protostars. In contrast, the sink 
particle in Run A is embedded in an accretion flow at all times. The ionizing radiation does 
blow away a significant fraction of gas from the sink particle, but it is not able to do so in 
all directions. There is always a region in the disk midplane around the sink particle where 
the gas is dense enough that evaporation by the ionizing radiation remains insufficient to 
stop the accretion flow. 



4. Bipolar Outflows 

In all runs, the H II regions are trapped in the disk plane but drive a bipolar outflow 
perpendicular to the disk. The highly variable rate of accretion onto protostars as they 
pass through dense filaments causes fast ionization and recombination of large parts of the 
interior of the perpendicular outfiow. The H ll regions around the massive protostars do 
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Fig. 5. — Disk fragmentation in the multiple sink simulation (Run B). The figure shows the 
gas density in slices in the midplane of the disk. The filaments in the disk form a succes- 
sively growing number of protostars. As the gas reservoir around the massive protostars is 
exhausted, the thermal pressure of the ionized gas creates a bubble around the star. This 
stops the accretion process. In the course of the simulation, the bubbles grow in size and 
finally merge. 
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Fig. 6. — Disk fragmentation in the single sink simulation (Run A). The images show the 
same region as Figure [5| This time the filaments do not form protostars but blobs of gas 
which are puffed up by the dynamical temperature fioor. The ionizing radiation is unable 
to create a bubble in the strong accretion fiow large enough to stop accretion. 
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not uniformly expand, but instead fluctuate rapidly in size, shape and flux. 

The diflFerences between the single and multiple fragment Runs A and Run B are also 
evident from the characteristics of the bipolar outflows. Figures [7| and |8] show density slices 
perpendicular to the disk plane for both simulations. Figure [7| shows the outflows driven by 
the most massive sink particle in Run B. The eflScient shielding by the fllaments and the 
motion of the sink particle hinder the thermal pressure of the ionized gas from driving a 
symmetric bipolar outflow at early times. The shielding of the ionizing radiation leads to a 
drop in the thermal pressure which drives and supports the outflow, and thus the outflow 
can be quenched again. At later times, the strong accretion flow onto the sink particle 
stops, and this allows the sink particle to drive a much larger outflow, which in the end has 
the form of a bubble. 

The situation for the sink particle in Run A is presented in Figure [8| The general 
properties of the outflow are very similar, but the big diflFerence is that here the accretion 
flow never stops. On the contrary. Figure |3] shows that the accretion rate increases towards 
the end of the simulation. This strong accretion flow fully absorbs the ionizing radiation, so 
the sink particle is unable to build up a bubble, even though it has more than three times 
the mass of the most massive star in Run B. 



Observed outflows ( Arce et aL||2007 ) from 0-type protostars with ages > 10^ yr tend to 



be much wider than those from low- or even intermediate-mass protostars. When resolved, 
these outflows show a disordered appearance that may be the outcome of multiple powering 
sources (including H ll regions) and non-steady driving. The pressure-driven outflows 
produced by H ll regions in our simulations reproduce these broad, disordered outflows 
well. Collimated outflows from less massive protostars without H ll regions are likely 



magnetically driven, though (see e.g., Banerjee & Pudritz 2006). A more detailed model 



linking these two types of outflows is out of the scope of this paper. 
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Fig. 7. — Bipolar outflows in the multiple sink simulation (Run B). The slices show the gas 
density in a plane through the most massive sink particle perpendicular to the disk. In the 
beginning of the simulation, the outflow is bipolar but not symmetric because the fllaments 
prevent the thermal pressure to drive simultaneously both lobes of the outflow. This eflFect 
fades when the accretion flow drops, and the sink particle can even blow an expanding bubble 
although it does not grow in mass anymore in the last three frames shown. The scales and 
times of the images are the same as in Figure [5) The slice plane is chosen such that it 
includes the sink particle of interest and is centered around it. 
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Fig. 8. — Bipolar outflows in the single sink simulation (Run A). The images show that the 
size of the outflow is not directly related to mass of the protostar. Instead, it is given by the 
time the thermal pressure has had to steadily drive it. If the thermal support is lost due to 
recombination, the outflow gets quenched again. This is apparent in the last frame, where 
the star is the most massive but there is no visible cavity around it. The scales and times of 
the images are the same as in Figure [6) The slice plane is chosen such that it includes the 
sink particle of interest and is centered around it. 
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5. Comparison to Observations 

We can directly compare our models with radio observations of free-free continuum, 
hydrogen recombination lines, and NH3(3,3) inversion lines by generating synthetic maps. 



5.1. Methods 



We generate radio continuum maps from our numerical data by integrating the 



radiative transfer equation in the Ray leigh- Jeans limit while neglecting scattering (Kraus 



1966; Gordon & Sorochenko 2002). To this end, we calculate for every cell in the domain 



the free-free absorption coefficient for atomic hydrogen 

-1.35 

a,, = 0.212 



-2.1 



cm 



(8) 



lcm-3y \1KJ VlHz 
where rie is the number density of free electrons, T the gas temperature and 1/ the frequency. 
The optical depth at distance r is then 



a^, ds. 



(9) 



Without scattering, the radiative transfer equation for the brightness temperature can be 
directly integrated and yields 



T(r.) = e-^'^ r e^nrl) d<. 

^0 



(10) 



If l^s is the solid angle subtended by the beam, the resulting flux density is 



A2 



111 



Following the algorithm described in Mac Low et al. (1991), we convolve the resulting image 



with the beam width and add some noise according to the telescope parameters. 

For a time-series comparison at high temporal resolution, we select a few time intervals 
in the multi-sink run (Run B) during which the continuum emission from a particular. 



-24- 



well-isolated, H ll region varies strongly. We then produced 2 cm continuum maps every 
^ 10 yr from the simulation output and measured the properties of the region, including 
the observational scale length i7, the flux F2cm, and the accretion rate of the protostar 
^sink that powers the H ll region. H is deflned as the equivalent diameter of a circle with 
the same area as the emission. F2cm is obtained by integrating the intensity of the maps 
over solid angle. Mgink is the accretion rate of the sink particle. 

We generate synthetic observations of line emission using our radiative transfer 



code MOLLIE (Keto 1990). This code generally uses the Accelerated Lambda iteration 



algorithm (Rybicki & Hummer 1991) although in this investigation we compute the level 
populations of NH3 assuming local thermodynamic equilibrium (LTE) and assume optically 
thin conditions for the H53a recombination line. The computation is done on 4 nested grids 
with spatial scales ranging from 80 to 640 AU and covering 0.2 pc. The maps of the NH3 
and H53Qf velocities show the flrst moment, in the case of NH3 integrating over the main 
hyperflne line of the (3,3) transition. 



5.2. Results 



Our simulated observations of radio continuum emission reproduce the morphologies 



reported in surveys of ultracompact H ll regions (Wood & Churchwell 1989; Kurtz et al. 



1994). Figured shows typical images from Run B, and the Online Material contains a 



movie of the entire run. The H ll regions in our model fluctuate rapidly between different 
shapes while accretion onto the protostar continues. When the gas reservoir around the two 
most massive stars is exhausted, their H ll regions merge into a compact H ll region similar 



to the extended envelopes often observed around ultracompact H ll regions (Kim & Koo 



2001; Kurtz et al. 1999). Our results suggest that the lifetime problem of ultracompact 



H II regions (Wood & Churchwell 19891) is only apparent. Since H ll regions embedded in 
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accretion flows are continuously fed, and since they flicker with variations in the flow rate, 
their size does not depend on their age until late in their lifetimes. 



The same mechanisms that cause the diflFerent accretion behavior (§ [3|) and bipolar 
outflows (§[4]) for Run A and Run B also lead to differences in the H ll region morphologies. 
One the one hand, the permanent accretion flow in Run A quenches the H ll region in 
the disk plane. On the other hand, the higher stellar mass in Run A leads to larger H ll 
regions perpendicular to the disk plane. Since the second eflFect dominates, the bottom 
line is that the H ll regions in Run A are generally larger than in Run B. Hence, multiple 
sink formation is not only necessary to model the gas fragmentation correctly, it is also 
crucial to reproduce the large relative fraction of spherical H ll regions observed in surveys 



(Wood & Churchwell 1989; Kurtz et al. 1994). A detailed investigation of the H ll region 



morpholgies, including a comprehensive statistical analysis for both simulations, will be 
presented in a follow-up paper. 

To compare more directly to observations of the time variability of H ll regions, we 
analyzed a few time intervals of interest at a time resolution of ^ 10 yr. We flnd that 
when the accretion rate to the star powering the H ll region has a large, sudden increase, 
the ionized region shrinks, and then slowly re-expands. This agrees with the contraction, 
changes in shape, or anisotropic expansion observed in radio continuum observations of 



ultracompact H ll regions over intervals of ^ 10 yr (Franco-Hernandez & Rodriguez 2004 



Rodriguez et al. 2007; Galvan-Madrid et al. 2008). Figure 10 shows that the sudden 



accretion of large amounts of material is accompanied by a fast decrease in the observed 
size and flux of the H ll region. In the left panel, the H ll region is initially relatively large, 
and accretion is almost shut off. A large (from to 6 x lO^^M© yr~^), sudden accretion 
event causes the H ll region to shrink and decrease in flux. The star at this moment has 
a mass of 19.8 Mq. In the right panel, the star has a larger mass (23.3 Mq), the H ll 
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Fig. 9. — H II region morphologies. This figure shows ultracompact H ll regions around 
massive protostars in Run B at different time steps and from different viewpoints. The 
cluster is assumed to be 2.65 kpc away, the full width at half maximum of the beam is 0^.^14 
and the noise level is 10~^Jy. This corresponds to typical VLA parameters at a wavelength 
of 2 cm. The protostellar mass of the central star which powers the H ll region is given in 
the images. The H ll region morphology is highly variable in time and shape, taking the 



form of any observed type (Wood & Churchwell 1989; Kurtz et al. 1994) during the cluster 
evolution. 



region is initially smaller, and the star is constantly accreting gas. The ionizing-photon 
flux appears to be able to ionize the infalling gas stably, until a peak in the accretion 
rate by a factor of three and the subsequent continuous accretion of gas makes the H ll 
region shrink and decrease in flux. The H ll region does not shrink immediately after the 
accretion peak because the increase is relatively mild and the geometry of the infalling gas 
permitted ionizing photons to escape in one direction. Our results show that observations 



of large, fast changes in ultracompact H ll regions (Franco-Hernandez & Rodriguez 2004 



Rodriguez et al. 2007; Galvan-Madrid et al. 2008) are controlled by the accretion process. 



Both the scale length and flux decrease at rates of 5-7 % per year in our model, agreeing 



well with observed fluctuations of 2~9 % per year (Franco-Hernandez & Rodriguez 2004 



Galvan-Madrid et al. 2008). Shortly after the minimum values are reached, the H ll regions 
re-expand, on timescales ^ 10^ yr. 

We can compare our models to a well-studied example of an ultracompact H ll region. 



W51e2 (Zhang et al. 1998; Keto & Klaassen 2008). Although we did not choose parameters 



speciflcally to model this region, the comparison is nevertheless revealing. In Figure 11 



we show simulated and observed maps of NH3(3,3) emission, 1.3 cm thermal continuum 
emission, and the H53a radio recombination line. The simulated maps were made at a time 
when a 20 Mq star has formed in Run B. 

In the simulated observations, the origin of the coordinates is referenced to the location 
of the 20 Mq protostar at the center of the accretion flow as marked. The accretion flow 
and disk are viewed edge-on. The spatial scale in the observations (lower panels) is 5100 
AU per arc sec assuming W51e2 is at a distance of 5.1 kpc ( |Xu et al.|[2009l ). The color 
bar on the right of each flgure shows the molecular velocities in km s~^. The velocities of 
the observations include the LSR velocity of W51e2, approximately 57 kms~^. The white 
contour levels are 50% through 90% of the peak brightness temperature of 71 K (upper left) 
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and 0.1, 0.2 and 0.3 Jy beam"-^ kms"-^ (lower left). The red contour levels are 30%, 70%, 
and 95% of the peak 1.3 cm continuum brightness temperature of 8902 K (upper left) and 
0.07, 0.14, and 0.22 Jy/beam (lower left). The white contours on the H53a observations 
(lower-right) show the 7 mm continuum emission at 2, 4, 10, 30, 50, 70, and 90% of the 



peak emission of 0.15 Jy beam ^. The molecular line observations are from Zhang et al. 



(1998). The H53a observations are from Keto & Klaassen (2008). Both observations have 



a beamsize of about V\ Coordinates are in the B1950 epoch. 

The brightest NH3(3,3) emission indicates the dense accretion disk surrounding the 
most massive star in the model, one of several within the larger-scale rotationally flattened 
flow. The disk shows the signature of rotation, a gradient from redshifted to blueshifted 
velocities across the star. A rotating accretion flow is identiflable in the observations, 
oriented from the SE (red velocities) to the NW (blue velocities) at a projection angle of 
135° east of north (counterclockwise). 

The 1.3 cm radio continuum traces the ionized gas, which in the model expands 
downward perpendicularly to the accretion disk, down the steepest density gradient. As 
a result, the simulated map shows the brightest radio continuum emission just off the 
mid-plane of the accretion disk, separated from the central star, rather than surrounding 
it spherically. In the observations of W51e2, continuum emission is indeed oflFset from the 
accretion disk traced in ammonia. The NH3(3,3) in front of the H ll region is seen in 
absorption and red-shifted by its inward flow toward the protostar. The density gradient in 
the ionized flow determines the apparent size of the H ll region. Therefore the accretion 
time scale determines the age of the H ll region rather than the much shorter sound-crossing 
time of the H ll region. 

Photoevaporation of the rotationally flattened molecular accretion flow supplies the 
ionized outflow. Therefore, the ionized gas rotates as it flows outward, tracing a spiral. 
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An observation that only partially resolves the spatial structure of the ionized flow sees a 
velocity gradient oriented in a direction between that of rotation and of outflow, as shown 



in the simulated observation. The observed H53Qf recombination line (Keto & Klaassen 



2008) in Figure 11 indeed shows a velocity gradient oriented between the directions of 



rotation and the outflow. 



6. Discussion 



6.1. Context 



In this work, we focus on the impact of feedback from radiative heating by both ionizing 
and non-ionizing radiation during massive star formation. Our work is complementary to 
recent simulations that neglect ionizing radiation but include other effects. 

For example, we use a ray-tracing algorithm to describe ionizing and non-ionizing 



radiation while Yorke fc Sonnhalter (2002); Krumholz et al. (2007b, 2009b); Price fc Bate 



(2009) rely on flux-limited diflFusion to follow the radiation fleld. Like us, these models 



include heating from accretion luminosity. They neglect heating by ionizing radiation, 
but do include radiation pressure and scattering effects. Radiation pressure in particular 
becomes important at the 10 AU scale. These models therefore simulate rather smaller 
scales than we do, starting from a ^ 100 core of only ^ 0.1 pc radius, as opposed to 
our more massive clump with 1000 Mq and a 1.6 pc radius. As a result, our calculations 
are better geared towards describing the cluster-forming scale as opposed to modeling in 



detail the mass growth of individual protostars through their inner accretion disks. Yorke 



& Sonnhalter (2002); Krumholz et al. (2007b, 2009b); Sigalotti et al. (2009) show that 



radiation pressure cannot stop accretion onto the massive protostar and is dynamically 
unimportant except for massive star clusters near the Galactic center or in starburst 
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galaxies (Krumholz & Matzner 2009). Krumholz et al. (2007b) initialized a microturbulent 



flow within their dense core, but Krumholz et al. (2009b) found that it made little difference 



by the time material had collapsed to the size of the disk. In addition, Price & Bate (2009) 
include magnetic fields and find that they provide some support for low-density gas, thus 



reducing the amount of matter available for accretion. Banerjee & Pudritz (2006), Machida 



et al. (2005), Hennebelle & Teyssier (2008), Hennebelle & Fromang (2008) and Hennebelle 



& Ciardi (2009) furthermore show that the presence of magnetic fields can dramatically 



alter the shape and size of accretion disks and strongly suppress the formation of close 
binary systems. These last results, may, however, be an artifact of the neglect of turbulent 



and ambipolar diffusion within the disks in these models. Dale & Bonnell (2008) and Wang 



et al. (2009), on the other hand, focused on mechanical feedback through winds and bipolar 



outfiows, without and with magnetic fields respectively, but neglect radiative heating 
entirely. They find that outfiows further limit accretion, but do not fundamentally change 
the picture of cluster formation in a gravitationally unstable accretion fiow. Ultimately, the 
quantitative characterization of the star formation efficiency will require treating the effects 
of radiation, magnetic fields, and outfiows simultaneously. 



6.2. Limitations 

As noted above, we have neglected a number of physical processes, including radiation 
pressure, scattering, magnetic fields, winds and bipolar outfiows, and turbulent initial 
conditions. In addition, we here discuss other approximations that we have made. 

Our prestellar model is designed such that the number of UV photons emitted by the 
protostar is overestimated for three reasons. First, we use a ZAMS model to set the ionizing 



luminosity, while the models of Hosokawa & Omukai (2009) suggest that the protostellar 



radii of massive stars can be larger than the ZAMS value, leading to a smaller ionizing 
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luminosity compared to the ZAMS value. Second, we neglect the aborption of UV photons 
by dust inside the H ll region. Thus, the sizes of our H ll regions represent upper limits to 
the real values. Third, the dense inner disk around each protostar is subsumed in the sink 
particle, allowing radiation to propagate along the disk midplane to the edge of the sink 
particle, where densities will be lower. 

However, as discussed in § [3| we find that, although we overestimate the strength of 
the ionizing radiation, it is unable to stop accretion onto a single massive star. Similarly, 
in § [5] we showed that the strength of the accretion flow restricts the expansion of the H II 



region, which solves the lifetime problem of ultracompact H ll regions (Wood & Churchwell 



1989). Both results are robust in the sense that the impact of the ionizing radiation on 
the accretion flow and the sizes of the H ll regions would even be smaller if we used a 
more detailed modeling of the protostellar evolution and environment, and included dust 
absorption. 

We also note that, although we nominally use a simplifled treatment of radiative 
losses in the optically thick regime, this will not have an impact on the dynamics of the 
simulated system, as we expect that the gas will become optically thick at densities of 



10 gem ^ (e.g. Larson 1969), which is two orders of magnitude above the threshold 



density for sink creation. 



7. Conclusions 

Our results suggest that the accretion flow onto massive stars cannot be stopped by 
ionizing radiation (§|3]). Instead, the growth of massive stars is halted by competition from 
lower-mass companions that form by gravitational fragmentation in the dense, rotationally 
flattened, accretion flow. This process, which we call fragmentation-induced starvation. 
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limits the maximum stellar mass of the highest-mass star in our simulations. 

Our models lead to a self-consistent picture of massive star formation where massive 
stars form through (unsteady) mass growth within a gravitationally unstable accretion flow. 
Radiation feedback from the accretion luminosity and the protostellar surface heats up 



the environment in which the flrst massive protostar forms (Krumholz et al. 2007b). This 
increases the local Jeans mass, so that subsequent fragments in the neigborhood of this star 
are relatively massive, naturally resulting in a cluster of massive stars. 



Unlike in the competitive accretion model (Bonnell & Bate 2006) where the most 
massive protostar continues to accrete from the common gas reservoir, in our model 
accretion is shut off by the subsequent fragmentation of the gas because the lower-mass 
companions intercept material that would otherwise accrete onto the central star, starving 
it of material. Fragmentation-induced starvation will only occur in regions with high enough 
density for many massive fragments to form despite accretion heating. In our models, the 
most massive star starts to form early and completes accretion while the whole cluster is 
still growing, contrary to standard competitive accretion models in which the most massive 
star grows over the whole cluster evolution timescale. On the other hand, competition for 
a common reservoir continues to deflne the top end of the initial mass function, unlike in 
models that rely on the clump mass spectrum to deflne the stellar initial mass function. 



Our models show (§ [5]) that the H ll regions formed around massive young stars flicker 
but do not grow systematically while heavy accretion from the protostellar core continues. 
This behavior seems able to explain the many puzzles raised by observations of ultracompact 
H II regions. The chaotic interaction between gravitationally unstable gas fllaments in the 
accretion disk and ionizing radiation from the young stars creates a great variety of different 
ultracompact H ll region morphologies. We flnd shell-like, core-halo, cometary, spherical 
and irregular morphologies in a single simulation, depending on the current flow fleld and 
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viewing angle. Our models show behavior consistent with observations of shrinking H ll 
regions with changes in size or flux of a few percent per year, as a consequence of denser 
material being accreted the star and shielding the ionizing radiation. Shielding by the disk 
also produces uncollimated bipolar outflows in our models (see § These results lead 
to an evolution scenario for H ll regions that is substantially different from the classical 
picture of monotonically expanding hot bubbles. The surprising dynamics of H II region 
flickering resolves the lifetime problem for ultracompact H ll regions, since the size of the 
H II region remains unrelated to the age of the protostar, as long as fast accretion continues. 
A detailed comparison of simulated observations of our model with the ultracompact H ll 
region W51e2 demonstrates that our model reproduces distinctive features of massive star 
forming regions such as H II regions expanding asymmetrically away from the protostar 
down the steepest density gradient, and out ward- twisting, spiraling, ionized flows. All 
these observations can be understood as consequences of ionization within a gravitationally 
unstable accretion flow. 
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Fig. 10. — Large amounts of molecular gas accreting onto an H ll region can cause a sudden 
decrease in its size and flux. This flgure shows the 2 cm continuum flux F2cm (in units of 
10 mJy), the characteristic size of the H ll region H (in units of 10 AU), and the rate of 
accretion to the star Mgink (sink particle, in units of lO^^M© yr~^) for two accretion events 
at high temporal resolution. 
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Fig. 11. — Comparison of line and continuum emission simulated from the model (upper 
panels) and actually observed from the W51e2 region (lower panels). The left panels show 
the NH3(3,3) line emission strength in white contours, the molecular line velocities as the 
background color, and the 1.3 cm free-free continuum from ionized gas in red contours. The 
right panels show the H53Qf recombination line velocities from the ionized gas. The simulated 
H53a observations (upper-right) are convolved to a spatial resolution of 750 AU (FWHM) to 
better match the relative resolution of the actual observations (lower-right). The molecular 



line observations are from Zhang et al. (1998). The H53a observations are from Keto & 



Klaassen (2008) 
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